Demonstration: Compare Density Models

This notebook demonstrates how to compare the GEODYN output for different density models

First, load the data using the reading functions.

%load_ext autoreload
%autoreload 2

import pandas as pd
import numpy as np

##################################
# INPUT PARAMETERS:
##################################
sat_file = 'starlette'
arc1 = '030914_2wk'
grav_id ='goco05s'
local_path = '/data/analysis/starlette_analysis/'
SAT_ID = 7501001

##################################
# PATH TO DENSITY MODEL RUN of Choice:
##################################
msis86_model = 'msis86'
path_to_msis86 = '/data/runs_geodyn/st/results/'+ grav_id+'_'+msis86_model+'/'

dtm87_model = 'dtm87'
path_to_dtm87 = '/data/runs_geodyn/st/results/'+ grav_id+'_'+dtm87_model+'/'

jaachia71_model = 'jaachia71'
path_to_jaachia71 = '/data/runs_geodyn/st/results/'+ grav_id+'_'+jaachia71_model+'/'


import sys
sys.path.insert(0, '/data/analysis/util_funcs/py_starlette/')
from a_ReadStarlette import ReadStarlette

(AdjustedParams_msis86,
 Trajectory_msis86,
 Density_msis86,
 Resids_msis86) = ReadStarlette(arc1,
                         sat_file,
                         grav_id,
                         local_path,
                         path_to_msis86)

(AdjustedParams_dtm87,
 Trajectory_dtm87,
 Density_dtm87,
 Resid_dtm87s) = ReadStarlette(arc1,
                         sat_file,
                         grav_id,
                         local_path,
                         path_to_dtm87)

(AdjustedParams_jaachia,
 Trajectory_jaachia,
 Density_jaachia,
 Resids_jaachia) = ReadStarlette(arc1,
                         sat_file,
                         grav_id,
                         local_path,
                         path_to_jaachia71)

The base file name for this arc is: st030914_2wk.goco05s

File exists: iieout
        /data/runs_geodyn/st/results/goco05s_msis86/IIEOUT/st030914_2wk.goco05s
File exists: ascii_xyz
        /data/runs_geodyn/st/results/goco05s_msis86/XYZ_TRAJ/st030914_2wk.goco05s
File exists: densityfil
        /data/runs_geodyn/st/results/goco05s_msis86/DENSITY/st030914_2wk.goco05s

 Loading data...

Parameter adjustment data loaded
Trajectory data loaded
Density data loaded
Residual data loaded
The base file name for this arc is: st030914_2wk.goco05s

File exists: iieout
        /data/runs_geodyn/st/results/goco05s_dtm87/IIEOUT/st030914_2wk.goco05s
File exists: ascii_xyz
        /data/runs_geodyn/st/results/goco05s_dtm87/XYZ_TRAJ/st030914_2wk.goco05s
File exists: densityfil
        /data/runs_geodyn/st/results/goco05s_dtm87/DENSITY/st030914_2wk.goco05s

 Loading data...

Parameter adjustment data loaded
Trajectory data loaded
Density data loaded
Residual data loaded
The base file name for this arc is: st030914_2wk.goco05s

File exists: iieout
        /data/runs_geodyn/st/results/goco05s_jaachia71/IIEOUT/st030914_2wk.goco05s
File exists: ascii_xyz
        /data/runs_geodyn/st/results/goco05s_jaachia71/XYZ_TRAJ/st030914_2wk.goco05s
File exists: densityfil
        /data/runs_geodyn/st/results/goco05s_jaachia71/DENSITY/st030914_2wk.goco05s

 Loading data...

Parameter adjustment data loaded
Trajectory data loaded
Density data loaded
Residual data loaded

Plot Drag Coefficient Comparison

labels = list(AdjustedParams_msis86[5]['0CD'].keys())
val_list = []
for i in AdjustedParams_msis86[5]['0CD'].keys():
    val_list.append(AdjustedParams_msis86[5]['0CD'][i]['CURRENT_VALUE'])

%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import rc
import matplotlib as mpl
from matplotlib import patches as mpatches
mpl.rcParams['lines.markersize'] = 10
plt.rcParams['axes.grid'] = True
plt.rcParams['grid.alpha'] = 1
plt.rcParams['grid.color'] = "#cccccc"
plt.rcParams.update({'font.size': 16})
rc('font',**{'family':'sans-serif','sans-serif':['Calibri']})
plt.rcParams['axes.titlesize']='large'
plt.rcParams['axes.titlepad']= 3
rc('text', usetex=False)
plt.rcParams["legend.loc"] = 'best'



color86 = 'tab:blue'
color71 = 'tab:red'
color87 = 'tab:green'

fig, ( ax1) = plt.subplots(1, figsize=(7,6), sharex=False)
for i in AdjustedParams_msis86.keys():
    ax1.set_title('Cd Output (Current Value) at each Iteration', y=1.1)
    ax1.plot(i, AdjustedParams_msis86[i]['0CD'][labels[0]]['CURRENT_VALUE'],'.', color = color86, label = 'MSIS 86 T01')
    ax1.plot(i, AdjustedParams_dtm87[i]['0CD'][labels[0]]['CURRENT_VALUE'],'x', color = color87, label = 'DTM T01')
    ax1.plot(i, AdjustedParams_jaachia[i]['0CD'][labels[0]]['CURRENT_VALUE'],'v', color = color71, label = 'Jaachia T01')


    ax1.set(ylabel=  'Cd [unitless]')
    ax1.set(xlabel=  'Iteration Number')
    msis_patch = mpatches.Patch(color=color86, label='MSIS')
    dtm_patch = mpatches.Patch(color=color87, label='DTM')
    jaachia_patch = mpatches.Patch(color=color71, label='Jaachia')

    fig.axes[0].legend( handles=[msis_patch, dtm_patch, jaachia_patch],bbox_to_anchor=(0, 1.0, 1.0, .08), ncol=3, mode="expand", borderaxespad=0.)

    fig.tight_layout( pad=1.0)


    for ax in fig.axes:
        plt.sca(ax)
        plt.xticks(rotation=45)

findfont: Font family ['sans-serif'] not found. Falling back to DejaVu Sans.
findfont: Font family ['sans-serif'] not found. Falling back to DejaVu Sans.

import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline
from plotly.subplots import make_subplots


last_iter = 5
which_stat = 'CURRENT_VALUE' # 2 is the current val

labels = list(AdjustedParams_msis86[5]['0CD'].keys())


fig = make_subplots(
    rows=1, cols=1,
    subplot_titles=("Time Dependent Drag Coefficient",
                    ))
val_list = []
for i in AdjustedParams_msis86[last_iter]['0CD'].keys():
    val_list.append(AdjustedParams_msis86[last_iter]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
                                 y=val_list,
                                  name= 'MSIS 86',
                               mode='markers',
                                    marker=dict(
                                    size=10,),
                                ), row=1, col=1,
                                   )

val_list = []
for i in AdjustedParams_dtm87[last_iter]['0CD'].keys():
    val_list.append(AdjustedParams_dtm87[last_iter]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
                                 y=val_list,
                                  name= 'DTM 87',
                               mode='markers',
                                    marker=dict(
                                    size=10,),
                                ), row=1, col=1,
                                   )

val_list = []
for i in AdjustedParams_jaachia[last_iter]['0CD'].keys():
    val_list.append(AdjustedParams_jaachia[last_iter]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
                                 y=val_list,
                                  name= 'Jaachia 71',
                               mode='markers',
                                    marker=dict(
                                    size=10,),
                                ), row=1, col=1,
                                   )

fig.update_yaxes( title=r"$C_D (Unitless)$",exponentformat= 'power',row=1, col=1)
fig.update_xaxes( title="Iteration Number", row=1, col=1)

# fig.update_yaxes( yaxis_title=r"$\frac{m}{s^2}$", xaxis_title="Iteration Number",exponentformat= 'power',row=2, col=1)

# fig.update_layout(
#     title="Along Track Accleration Amplitude vs Iteration",    )
# fig.update_layout(
#     yaxis_title=r"$\frac{m}{s^2}$",
#     xaxis_title="Iteration Number" )


iplot(fig)

Density along the orbit:

data_nums = 1000

import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline

fig = go.Figure(data=[go.Scatter(x=Density_msis86['Date'][:data_nums],
                                 y=Density_msis86['rho (kg/m**3)'][:data_nums],
                                    name='MSIS 86',
                                     mode='markers',
                                    marker=dict(
                                    size=2,
                                    ),
                                   ),
                                   ],
                                   )


fig.add_trace(go.Scatter(x=Density_dtm87['Date'][:data_nums],
                                 y=Density_dtm87['rho (kg/m**3)'][:data_nums],
                                  name= 'DTM 87',
                               mode='markers',
                                    marker=dict(
                                    size=2,),
                                ),
                                   )
fig.add_trace(go.Scatter(x=Density_jaachia['Date'][:data_nums],
                         y=Density_jaachia['rho (kg/m**3)'][:data_nums],
                         name= 'Jaachia 71',
                         mode='markers',
                         marker=dict(
                         size=2,),
                         ),
                         )

fig.update_layout(
    title="Density along Starlette Orbit (for 4 minutes of first day)",
    yaxis_title=r'$\frac{kg}{m^3}$',
    xaxis_title="Date",
    )
fig.update_yaxes(type="log", exponentformat= 'power',)

iplot(fig)

data_nums = 15
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline

fig = go.Figure(data=[go.Scatter(x=Density_msis86['Date'][::data_nums],
                                 y=Density_msis86['rho (kg/m**3)'][::data_nums],
                                 name = 'MSIS 86',
                                 mode='markers',
                                    marker=dict(
                                    size=2,
                                    ),
                                   ),
                                   ],
                                   )

fig.add_trace(go.Scatter(x=Density_dtm87['Date'][::data_nums],
                                 y=Density_dtm87['rho (kg/m**3)'][::data_nums],
                                  name= 'DTM 87',
                               mode='markers',
                                    marker=dict(
                                    size=2,),
                                ),
                                   )
fig.add_trace(go.Scatter(x=Density_jaachia['Date'][::data_nums],
                         y=Density_jaachia['rho (kg/m**3)'][::data_nums],
                         name= 'Jaachia 71',
                         mode='markers',
                         marker=dict(
                         size=2,),
                         ),
                         )
fig.update_layout(
    title="Density along Starlette Orbit (Every 15th datapoint)",
    yaxis_title=r'$\frac{kg}{m^3}$',
    xaxis_title="Date",
    )
fig.update_yaxes(type="log", exponentformat= 'power',)

iplot(fig)

import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline

fig = go.Figure(data=[go.Scatter(x=pd.to_datetime(Resids_msis86['Date']),
                                 y=Resids_msis86['Residual'].values.astype(float),
                                 name ='MSIS 86' ,
                                 mode='markers',
                                    marker=dict(
                                    size=3,
                                    ),
                                   ),
                                   ],
                                   )


fig.add_trace(go.Scatter(x=pd.to_datetime(Resid_dtm87s['Date']),
                                 y=Resid_dtm87s['Residual'].values.astype(float),
                                  name= 'DTM 87',
                               mode='markers',
                                    marker=dict(
                                    size=3,),
                                ),
                                   )
fig.add_trace(go.Scatter(x=pd.to_datetime(Resids_jaachia['Date']),
                         y=Resids_jaachia['Residual'].values.astype(float),
                         name= 'Jaachia 71',
                         mode='markers',
                         marker=dict(
                         size=3,),
                         ),
                         )

fig.update_layout(
    title="Observation Residuals from Final Iteration",
    yaxis_title="Residuals",
    xaxis_title="Date",
    )

iplot(fig)